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Abstract 



A time-domain numerical modeling of Biot poroelastic waves is presented. 
The viscous dissipation occurring in the pores is described using the dynamic 
permeability model developed by Johnson-Koplik-Dashen (JKD). Some of 
the coefficients in the Biot-JKD model are proportional to the square root 
of the frequency: in the time-domain, these coefficients introduce order 1/2 
shifted fractional derivatives involving a convolution product. Based on a dif- 
fusive representation, the convolution kernel is replaced by a finite number 
of memory variables that satisfy local-in-time ordinary differential equations. 
Thanks to the dispersion relation, the coefficients in the diffusive representa- 
tion are obtained by performing an optimization procedure in the frequency 
range of interest. A splitting strategy is then applied numerically: the prop- 
agative part of Biot-JKD equations is discretized using a fourth-order ADER 
scheme on a Cartesian grid, whereas the diffusive part is solved exactly. Com- 
parisons with analytical solutions show the efficiency and the accuracy of this 
approach. 
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1. Introduction 



Porous media consist of a solid matrix within which fluids can circulate 
freely. The propagation of waves in these media has many crucial implica- 
tions in applied mechanics, in situations where materials such as industrial 
foams, spongious bones j32[ and petroleum rocks 0] have to be characterized, 
for example. The poroelastic model originally developed by Biot in 1956 [l| 
includes two classical waves (one "fast" compressional wave and one shear 
wave), in addition to a second "slow" compressional wave, which is highly de- 
pendent on the saturating fluid. This slow wave was observed experimentally 



in 1981 30], thus confirming the validity of Biot's theory. 

Two frequency regimes have to be distinguished when dealing with poroe- 
lastic waves. One of the main problems is how to model the viscous efforts 
generated bv the mechanical perturbations involved. In the low-frequency 
range (LF) the viscous boundary layer that develops in the fluid is large 
in comparison with the diameter of the pores, and the viscous efforts are pro- 
portional to the relative velocity of the motion between the fluid and solid 
components. In the high-frequency range (HF), modeling the dissipation is 
a more delicate task: Biot first presented an expression for particular pore 
geometries In 1987, Johnson-Koplik-Dashen (JKD) [l9j published a gen- 
eral expression for the dissipation in the case of random pores. The viscous 
efforts depend in this model on the square root of the frequency of the per- 
turbation. When writing the evolution equations in the time domain, time 
fractional derivatives are introduced, which involves convolution products 
with singular kernels 25)]. Analytical solutions have been derived in simple 



academic geometries and homogeneous media [13 



Many numerical methods have been developed in the LF regime: see |5| 
and the introduction to Q for general reviews. In the HF regime, the frac- 
tional derivatives greatly complicate the numerical modeling of the Biot- JKD 
equations. The past values of the solution are indeed required in order to 
evaluate these convolution products, which means that the time evolution of 
the solution must be stored. This of course greatly increases the memory re- 
quirements and makes large-scale simulations impossible. To our knowledge, 
only two approaches to this problem have been proposed so far in the liter- 
ature. The first approach consisted in discretizing the convolution products 



26j . and the second one was based on the use of a diffusive representation 
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of the fractional derivative [24], |34|. In the latter approach, the convolution 
product is replaced by a continuum of diffusive variables - or memory vari- 
ables - satisfying local differential equations [17]. This continuum is then 
discretized using appropriate quadrature formulas, resulting in the Biot-DA 
(diffusive approximation) model. 

However, the diffusive approximation proposed in 24j has three major 
drawbacks. First, the quadrature formulas make the convergence towards 
the original fractional operator very slow. Secondly, in the case of small fre- 
quencies, the Biot-DA model does not converge towards the Biot-LF model. 
Lastly, the number of memory variables required is not specified. The aim 
of the present study is therefore to develop a new diffusive approximation 
method in which these drawbacks do not arise. Since it is proposed here to 
focus on the discretization of the fractional derivatives, we will deal only with 
the 1-D equations of evolution in homogeneous media, so that the shear wave 
will not be considered. However, the strategy proposed here can be extended 
quite straightforwardly to 2D and 3D geometries, as discussed below. 

This paper is organized as follows. The original Biot-JKD model is briefly 
outlined in section [2] and the principles underlying the diffusive representa- 
tion of fractional derivatives are described. The decrease of energy and the 
dispersion analysis are addressed. In section [3j the method used to discretize 
the diffusive model is presented: the diffusive approximation thus obtained is 
easily treatable by computers. Following a similar approach than in viscoelas- 



ticity [15] , the coefficients of the model are determined using an optimization 
procedure in the frequency range of interest, giving an optimum number of 
additional computational arrays. The numerical modeling is addressed in 
section HI where the equations of evolution are split into two parts: a prop- 
agative part, which is discretized using a fourth-order scheme for hyperbolic 
equations, and a diffusive part, which is solved exactly. Some numerical 
experiments performed with realistic values of the physical parameters are 
presented in section |5j In section [6j a conclusion is drawn and some futures 
lines of research are given. 



2. Physical modeling 

2.1. Biot model 

The Biot model describes the propagation of mechanical waves in a macro- 
scopic porous medium consisting of a solid matrix saturated with a fluid 
circulating freely through the pores [H, S 0] . It is assumed that 
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• the wavelengths are large in comparison with the diameter of the pores; 



• the amplitude of the perturbations is small; 

• the elastic and isotropic matrix is completely saturated with a single 
fluid phase; 

• the thermo-mechanical effects are neglected. 

This model involves 10 physical parameters: the density pj and the dynamic 
viscosity 77 of the fluid; the density p s and the shear modulus p of the elastic 
skeleton; the porosity < < 1, the tortuosity a > 1, the absolute perme- 
ability k, the Lame coefficient A/ and the two Biot's coefficients and m of 
the saturated matrix. The following notations are introduced 

Pw = ^Pf, P = <j> Pi + (1 - 4>) Ps: X = PPw-p}>0, 

A = A/-m/3 2 , C = A + 2/i>0. 

Taking u s and Uf to denote the solid and fluid displacements, the unknowns 
in ID are the elastic velocity v s = ^f, the filtration velocity w = ^ff = 
(j) £ (uf — u s ), the elastic stress a, and the acoustic pressure p. The consti- 
tutive laws are 

a = (\ f + 2n)e-mP£, (2a) 
p = TO (_/3 £ + £) ) (2b) 
where e = is the strain and £ = is the rate of fluid change. On the 

ox ^ ox 

other hand, the conservation of momentum yields 

dv s dw da . . 

P ^T + P ^ = ^ (3a) 

dv s dw r] dp 

Ps ^r +Pw ^ + - K F * W = -^ (3b) 



where * is the convolution product in time. The equation (I3bl) is a generalized 
Darcy law. The quantity F * w denotes the viscous dissipation induced by 
the relative motion between the fluid and the elastic skeleton. 
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2.2. High frequency dissipation: the JKD model 

The frontier between the low-frequency range (LF) and the high-frequency 
range (HF) is reached when the viscous efforts and the inertial effects are 
similar. The transition frequency is given by [H, \^ 

h= -H^ = ^. (4) 

In LF, the flow in the pores is of the Poiseuille type, and dissipation efforts 
in ( |3b| are given by 

F(t) = 5(t) F(t) * w(x, t) = w(x, t), (5) 

where 5 is the Dirac distribution. In HF, the width of the viscous boundary- 
layer is small in comparison with the size of the pores, and modeling the 
dissipation process is a more complex task. Here we adopt the widely-used 
model proposed by Johnson-Koplik-Dashen (JKD) in 1987, which is valid 
for random networks of pores with constant radii [19j . The only additional 
parameter is the viscous characteristic length A. We take 

0A 2 ' P 4o 2 « 2 p/' 

where P is the Pride number (typically P ~ 1/2). Based on the Fourier 
transform in time, F(u) = j R F{t)e~ luJt dt, the frequency correction given by 
the JKD model can be written 

/ A 2 2 \ 1/2 

- . / Aa z K z p f ^ 



r)A 2 < 



■,2 



l + i P^) 1/ \ W 



1 



(VL + iuj) 1 ' 2 . 



This correction is the sim ples t function satisfying the LF and HF limits of 
the dynamic permeability 19]. Therefore, the term F(t) *w(x,t) involved in 
dSbJ) is 

F(t)*w(x,t) =T~ X ( (£1 + iu) 1/2 w(x, u 



n ' (8) 

(D + Q) 1/2 w{x,t). 
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The operator D 1 / 2 is a shifted order 1/2 time fractional derivative, generaliz- 
ing the usual derivative characterized by = (icow(uj)). The notation 
(D + Q) 1 / 2 accounts for the shift in ©. 

2.3. The Biot-JKD equations of evolution 

Based on 02]), (E]) and ([H]), the Biot-JKD equations can be written 



( dv x 



dw da 
8t + Pf l)t = ~dx~' 



Pf 



dv. 



+ Pv 



dw i] 1 



+ 



dt ^ dt K^U 

a = (Xf + 2 fi)e — + s 

lp = m (-(3 e + £). 



(D + O) 



V2, 



w = — 



dp 



(9a) 

(9b) 

(9c) 
(9d) 



We rearrange this system by separating ^ and in (pal) and (l9b| and 
using the definitions of e and £. Taking 
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V P 1 



one obtains the following system of equations of evolution 



at x dx x dx p 



(10) 



dw pf da p dp 
dt x dx x dx 



^-(A/ + 2/,)^ 



m I 3 1T~ = f° 
dx 



dp n dv s dw 

— + m(3—— + m— = f p . 

dt dx dx 



(11a) 
(lib) 
(11c) 
(lid) 



Terms f Vs , f w , f a and f p have been added to the equations of evolution to 
simulate sources. 



6 



2-4- The diffusive representation 
Taking 

w 

D n w{x,t) = — + Qw, (12) 
the shifted fractional derivative (jSJ) can be written [10] 

i ft -fi(t-r) 

p + fl) 1 / 2 «>(ar J t) = -= / ; D n w(x, r) dr. (13) 



The operator (D + Vt) 1 ^ 2 is not local in time and involves the entire time 
history of w. As we will see in section [3j a different way of writing this 
derivative is more convenient for numerical evaluation. Based on Euler's V 
function, the diffusive representation of the totally monotone function -t= 

@, 03, E3, is 

ii roo i 

- et d9. (14) 



Substituting (fl4"l) into ({TBI gives 

1 ft poo 1 

(D + fi) 1/2 w(x,t) =- / / -^e- ^e- a ^D n w(x,T)d0dr, 

k Jo Jo Vd 

1 f°° 1 

= - / -j=^(x,e,t)de, 

n Jo V0 

(15) 

where the diffusive variable is defined as 

1>(x,e,t)= [ e- ie+m - T) D n w(x,r)dT. (16) 
Jo 

For the sake of clarity, the dependence on fl and w is omitted in ip. From 
(ITS]) , it follows that the diffusive variable ^ satisfies the ordinary differential 
equation 

§* = -<« + + (17) 
ip(x,e,o) = o. 

The diffusive representation therefore transforms a non-local problem ([TBI 
into a continuum of local problems f|T7|) . It should be emphasized at this 
point that no approximations have been made up to now. The computational 
advantages of the diffusive representation will be seen in sections [3] and [51 
where the discretization of (|15|) and (|17|) will yield a tractable formulation. 
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2.5. Energy of Biot-JKD 

Now, we express the energy of the Biot-JKD model [9]). This result gen- 



eralizes the analysis performed in the LF range in [12 
Proposition 1. Let 

E = E\ + E 2 + i?3, 

with 

Ei = - (pv 2 + p w w 2 + 2p f vw) dx, 
2 Jm. 

W/l, - x2 1 



^2 = 2 J R \C {a + Pp) + m P > dX ' (18) 
E 3 = - f f l-^= n 1 - (w-xpfdBdx. 



Then E is an energy which satisfies 
dE f f i] 1 1 1 



(VLw 2 + (6 + tt)ip 2 )d6dx < 0. (19) 



Proposition [T] is proven in appendix 1. It calls for the following comments: 

• the Biot-JKD model is well-posed; 

• when the viscosity of the saturating fluid is neglected (77 = 0), the 
energy of the system is conserved; 

• the terms in f fl8|) have a clearly physical significance: E\ is the kinetic 
energy, and E 2 is the potential energy. The term E 3 corresponds to the 
kinetic energy resulting from the relative motion of the fluid inside the 
pores. 

2.6. Dispersion analysis 

Injecting a mode e. l ( ut ~ kx > in ( ITT1) gives the dispersion relation between 
the angular frequency u and the wavenumber k. Taking 

' £> 4 =m(\ + 2p), 
D 2 (co) = - ((A/ + 2 y) p w + m{p - 2 p f $)) u 2 + i u - F(u) (A/ + 2p), 

K 

D (u) =xw 4 -w 3 -pfM- 

K 

(20) 
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the dispersion relation takes the form 

D e (k, to) = Dik 4 + D 2 (lo) k 2 + D (co) = 0. (21) 

Expressions (I20l) - (I2T1) are valid in the case of both the Biot-LF and Biot-JKD 
models with the frequency correction defined by 

{F LF (co) = 1 Biot-LF, (22a) 

1 
F JKD (u) = —(Q + i co) 1/2 Biot-JKD. (22b) 

The solutions k p f and k ps of ( l2Tj) give the phase velocities c p f = uj /^R.e{k p f) 
of the fast wave and c ps = uo/^te{k ps ) of the slow wave, with < c ps < 
c p f. The attenuations a p f = —Qm(k p f) and a ps = -3m(fc ps ) can also be 
deduced. Both the phase velocities and the attenuations of Biot-LF and Biot- 
JKD are strictly increasing functions of the frequency. The high frequency 
limits of fast and slow phase velocities, c^- and c^, which are obtained by 
diagonalizing the left-hand side of system (fTTj) . satisfy the relation 

X c A - ((\f + 2 f i)p w + m(p-2p f P)) c 2 + m(\ + 2fi) = 0. (23) 

Figure [1] shows the dispersion curves corresponding to the Biot-LF and Biot- 
JKD models. The physical parameters are those used in the numerical ex- 
periments presented in section Note that the scales are radically different 
in the case of fast and slow waves. The following properties can be observed: 

• when / < f c , the Biot-JKD and Biot-LF dispersion curves are very 
similar as might be expected, since lim Fjkd(u) = 1; 

• the fast wave is almost not affected by the frequency correction F(u) 
while the slow wave is greatly affected; 

• when / f c , the slow wave is almost static. When / > f c , the slow 
wave propagates but is greatly attenuated. 

3. The Biot-DA (diffusive approximation) model 

The aim of this section is to approximate the Biot-JKD model, using a 
numerically tractable approach. 
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phase velocity of the fast wave 



phase velocity of the slow wave 
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3.1. The Biot-DA first-order system 

Using a quadrature formula on N points, with weights ag and abscissa 
0£ > 0, the diffusive representation (|T5|) can be approximated by 

1 f°° 1 

(D + Q) 1/2 w(x,t) =- / -=ip(x,t,0)d6 



^2a e ip(x,t,9t), (24) 
e=i 

N 

^2a e ipe(x,t). 



N 



1=1 

From (ITTj) . the iV diffusive variables ^ satisfy the ordinary differential equa- 
tions 

?£ = _ {0l + mi + DtilOt 
^(x,0) =0. 

The fractional derivatives are replaced by their diffusive approximation fl24l) 
in the JKD model ( fTTI) . Upon adding the equations ( l25l) and performing 
some straightforward operations, the Biot-DA system is written as a first- 
order system in time and space 

a N 

dv s 



Pw da p f dp Pf i , t 

■), » "F = TL^ + /«- 



A? 



<9w pf <9cr p dp v-^ , , 

Ot Y OX Y OX *— ' 

A A £=1 

_-(A / + 2„)^- m/3 - = /„, 

<9p dv s dw 
— + m/3 — + m — = f p , 
ot ox ox 

^ + P ^^ + ^-^-if:^e-(e j + n)i; J + f w , j ^i,...,N. 

A A 

(26) 

Taking the vector of unknowns 

U = (v s ,w,a,p, ipi, . . . ,?p N ) T (27) 



11 



and the source vector 



^ {fv s i fwi fai fpi fwi ■ ■ ■ i fw) 



the system (I26p can be written 



dXJ k d\J 
+ A 



dt d x 
where A is the (N + 4) 2 propagation matrix 



SU + F, 
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and S is the (N + 4) 2 dissipation matrix 



(° 
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(31) 

The size of the system increases linearly with the number iV of diffusive 
variables. 



3.2. Properties 

Four properties of system (!29l) are specified: 

• the eigenvalues of A (!30l) are real: with multiplicity N, ±c^ and ±c^, 
where the latter satisfies (1231) . The system (|2T?|) is therefore hyperbolic; 

• since the eigenvalues and eigenvectors do not depend on the diffusive 
coefficients, they are the same in both the Biot-DA and Biot-LF or 
Biot-JKD models. This is not so in the case of the method presented 
in 26(|, where the propagation matrix is modified to account for the 



fractional derivative; 

the dispersion analysis is obtained in the case of the Biot-DA model by 
replacing F by 



n , • N 
F AD {u) = f=- 2^ 



(If; 



^Ot + n + iu 



(32) 



in equations fl20j> - fl2l|> : 
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• in line with proposition [TJ an energy analysis of ( 12 6 p is performed. 
Proposition 2. Let 



with 



E = Ei + E 2 + E 3 , 

E t — ^ I (pv 2 + p w w 2 + 2pfvw) dx, 
* Jr 

1 f N 1 ^ 
Es — - 7 f= ^ 7T (w — ?/v) 2 <ir. 



T/ien satisfies 



dE f 7] 1 



Since the proof is very similar in this case, it will not be repeated. The 
terms E\ and E% are the same in both the Biot-DA and Biot-JKD models, 
whereas E3 and the time evolution of E differ; in Biot-DA, the sign depends 
on the coefficients introduced into the diffusive approximation. The abscissas 
Qi of the quadrature formula are positive, but no sign criterion is given a 
priori for the weights ag. E therefore cannot be said to be a decreasing 
energy, except in the obvious case where all the an are positive. 

3.3. Determining the Biot-DA parameters 



The a£ and 6e in ( 124)) now have to be determined. In [24|, the authors used 
a general Laguerre quadrature formulas. We have tried using this approach, 
but it gave poor results. Very large numbers of diffusive variables were re- 
quired to approximate the Biot-JKD model accurately, resulting in a huge 
computational cost. In addition, the Biot-DA model based on Laguerre func- 
tions does not converge by construction towards Biot-LF when the frequency 
tends towards 0, which is neither satisfactory nor physically realistic. Lastly, 
the involved coefficients do not depend on the physical factors (parameters, 
source) involved, which partly explains the above two weaknesses. 
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A different method of determining the 2N coefficients and 9# in the 
diffusive approximation (1251) is therefore used, in order to approach F JKD (u) 
(!22|) by Fad(u) (|32|) in a given frequency range of interest. Let Q{u) be the 
optimized quantity and Q re f(u) be the desired quantity 



Q{u) 



Fad(co) 
Fjkd{u) 



^ e e + n + i oj = ^ aMujh (35a) 



k Qref(0j) = 1- 



(35b) 



We implement a linear optimization procedure [ll|, [16|, [23j in order to min- 
imize the distance between Q(u) and Q re f(u) in the interval [uj m i n ,u max ] 
centered on ojq = 2 n f , where fo is the central frequency of the source. The 
abscissas 0£ are fixed and distributed linearly on a logarithmic scale 
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UJr, 



e-i 

N-l 



N. 



(36) 



The weights ae are obtained by solving the system 

N 



^2a e q E (u k ) = I, k = 1, K, 



(37) 



where the Uk are also distributed linearly on a logarithmic scale of K points 



Uk — UJ n 



0O r . 



OJr, 



fc-1 
K-l 



k = 1,...,K. 



(38) 



Since the qe(oo) are complex functions, optimization is performed simultane- 
ously on the real and imaginary parts 



( N 



^2a e Re(q e (u k )) =1, 

N 

^a^Im(g £ (o; fc )) =0, k = l,...,K. 



(39) 



1=1 



A square system is obtained when 2K = N, whereas 2 K > N yields an 
overdetermined system, which can be solved by writing normal equations 
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2K = N 



K = N 
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14|. For practical purposes, we use u min = uj /10 and co max = 10oj o , as in 
23|. 

Figure |2] illustrates the influence of N and K on the accuracy of the opti- 
mization procedure. As can be observed in this figure, the errors are smaller 
with the overdeterminated system (K = N, 2N, 3N) than with the square 
one. However, increasing the size of the system does not really improve the 
accuracy. In what follows, we will therefore always use the values K = N. 
The influence of the number of diffusive variables on the physical properties 
of the system is presented in figure [3j We focus here on the slow wave, since it 
is more sensitive to the frequency correction. As was to be expected, the ac- 
curacy of the approximation of the Biot-JKD phase velocity and attenuation 
given by the Biot-DA model increases with N. 




Figure 3: Phase velocity c ps (left) and attenuation a ps (right) of the slow wave obtained 
with the Biot-DA model, in terms of the number of diffusive variables. 

To determine N in terms of the required accuracy, e m = \\Q(co) — 1||l 2 is 
measured in the frequency range of interest [fo/10, 10/o]. This norm amounts 
to the relative error between Fad{w) and Fjkz>(uj). With N < 20, this error 
is proportional to iV -1 , as can be seen from figure H]- (a). At larger values of N, 
the system is poorly conditioned and the order of convergence deteriorates 
(not shown here); in practice, this is not penalizing, however, since large 
values of N are of no use. An example of the parametric determination of 
TV in terms of the frequency range and the desired accuracy is also given 
in figure IH-(b). In the following numerical tests, N = 6 variables are used, 
giving the modeling error e m ~ 5%. 
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Figure 4: Determining the number of diffusive variables N. Left: relative error e m in 
terms of the number of N; the dashed line is proportional to N . Right: required value 
of N in terms of /o// c and the required accuracy e m (b). 



Lastly, the sign of weights ag was examined in a large number of config- 
urations. In each case, some negative values were obtained with the linear 
optimization process ( |39l) . As stated in proposition El the well-possedness of 
Biot-DA could not therefore be proved. A nonlinear optimization procedure 
with a positivity constraint was then applied [29j . but almost all the at ob- 
tained were equal to zero. In the numerical experiments, the negativity of 
some ai has never raised any problems. This question is addressed in detail 
at the end of section I5.2I 



4. Numerical modeling 

4-1- Splitting 

In order to integrate the Biot-DA system ( 129|) . a uniform grid is intro- 
duced, with mesh size Ax and time step At. The approximation of the 
exact solution U(xj = j Ax,t n = n At) is denoted by U™. If an unsplit 
integration of ( 1291) is performed, Von-Neumann analysis typically yields the 
stability condition 

At<min|T— , -Ar | , (40) 
^ dfi R(S) J' V ' 

where R(S) is the spectral radius of S, and T > depends on the numerical 
scheme. We have no theoretical estimate of R(S), but numerical studies have 
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shown that this value is similar to that of the spectral radius in LF: - £ , which 

r-i 1 K x' 

can be very large [7j. The time step can therefore be highly penalized in this 
case flUD- 

A more efficient strategy is adopted here, which consists in splitting the 
original system ( 129]) into a propagative part and a diffusive part ( 142]) 



+ A- 



< 



dt 

dV 
~dt 



dx 



-SU. 



(41) 



(42) 



For the sake of simplicity, the source term F has been omitted here. The 
discrete operators associated with steps (f4T|) and ( l4"2l are denoted by H a 
and Hf,, respectively The second-order Strang splitting [21] is then used to 
integrate fl29|) between t n ant £ n +i, giving the time- marching 



u[ 2) 



ur 1 



H b (f)U?, 
H a (At) Uf, 



(43) 



H 6 (^)U^ 



(2) 



The discrete operator H a associated with the propagative part fj4Tl) is an 
ADER 4 (Arbitrary DERivatives) scheme [3l|. This scheme is fourth-order 
accurate in space and time, is dispersive of order 4 and dissipative of order 
6, and has a stability limit T = 1. On Cartesian grids, ADER 4 amounts to 
a fourth-order Lax-Wendroff scheme, and can be written 



H„(Ai)U 



E 



(1) 



7n 



771=1 



A 



£ c,uf>„ 

s=-2 



At 



Ax 



(44) 



where the coefficients 7 TO)S are given in table [TJ 

Since the physical parameters do not vary with time, the diffusive part 
( l4"2l) can be solved exactly. This gives 



At 



(45) 
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lm,s 


m = 1 


m = 2 


m = 3 


m = 4 


S = 


-2 


1/12 


1/24 


-1/12 


-1/24 


S = 


-1 


-2/3 


-2/3 


1/6 


1/6 


s = 








5/4 





-1/4 


s = 


+1 


2/3 


-2/3 


-1/6 


1/6 


s = 


+2 


-1/12 


1/24 


1/12 


-1/24 



Table 1: Coefficients of the ADER 4 scheme. 



The matrix e 2 s is computed numerically using the [r / q] Pade approxi- 



mation in the "scaling and squaring method" [28], which is given by the 
expression 



R 



rq 



N, 



rq 



D 



rq 



2 



2 



E 

fc=0 

g 

E 



iV, 



At 



S) 

(r + q — k) ! r ! 



fc=0 



(r + g) ! k ! (r — k) ! 

(r + g — fc) ! g ! 
(r + q) \ k \ (r — k)\ 



At 



2 



(46) 



In the following numerical experiments, the parameters r = q = 6 are used. 
It remains to verify that the numerical integration of the diffusive step 
15]) is unconditionally stable. This is achieved as follows. 



Proposition 3. The diffusive part of the splitting U]Q is well-posed whatever 
the weights in the diffusive approximation (24\ ). 

Proposition [3] is proven in appendix 2. It follows that the solution of 
system (IB.lj) is bounded and that the eigenvalues of — ^ S are then in the 
left half space. As a consequence, the R qq Pade approximation is always 
stable 28) . The full algorithm (j4"3"|) is therefore stable under the optimum 
stability condition 

At<T^, (47) 

c pf 
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which is always independent of the Biot-DA model coefficients. 

5. Numerical experiments 

5.1. General configuration 



Saturating fluid 


p f (kg/m 6 ) 


1040 




Ti i Pft t; 1 


1 5 1 o -3 


drain 


n„ (keVin 3 ") 
Ps V 6/ / 


2650 




ji (Pa) 


2.93 10 9 


Matrix 





3.35 lO" 1 




a 


2 




k (m 2 ) 


io- n 




A/ (Pa) 


6.14 10 9 




m (Pa) 


6.49 10 9 




/3 


9.56 lO" 1 




A (m) 


2.19 10~ 5 


Phase velocities 


<%f ( m / s ) 


2384.17 




c~ (m/s) 


758.95 




„oo /„oo 
c pf/ c ps 


3.14 




fc (Hz) 


3844.97 



Table 2: Physical parameters used in numerical experiments. 

The physical parameters used in all the numerical experiments, which 
are given in table |5J correspond to sandstone saturated with water. The 
unbounded medium is excited by a point source f a = g(t) h(x), with h(x) = 
8(x) in equation flllcp . The time-dependent evolution of the source, g(t) in 
(lllcp . is a C 6 combination of truncated sinusoids 

{21 63 1 
sin (con t) sin (2 u 1) H sin (4 uj t) sin (8 w 1) if < t < 
v ' 32 v 1 768 V ' 512 v u ' ~ ~ 
otherwise, 

(48) 

with a central frequency /o = = 200 kHz. Adopting the high-frequency 
regime is therefore completely justified since /o — 50 x f c . Figure [5] shows 
the time-dependent evolution and spectrum of the source. 
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10 4 
f(Hz) 



Figure 5: Time-dependent evolution (left) and spectrum (right) of the source. 



The computational domain [—0.04, 0.04] m is discretized with N x grid 
points, and the time step is deduced from (1471) . taking T = 0.9. No special 
care is applied to simulate outgoing waves (with PML, for instance), since 
the simulations are stopped before the waves have reached the edges of the 
computational domain. The numerical experiments are performed on an Intel 
Core i7 processor at 2.80 GHz. 

Exact solutions of time-domain Biot-JKD equations have been derived 



in the literature [13], but not for Biot-DA. Therefore, we compute reference 
solutions of both Biot-JKD and Biot-DA thanks to standard tools of Fourier 
analysis: the Green functions of (flTT) or (126]) are determined in the harmonic 
regime. Then, the Cauchy residue theorem and numerical inverse Fourier 
transforms (Nf = 9.6 10 5 modes and a frequency step A / ~ 13 Hz) yield the 
semi-analytical solutions. 

5.2. Test 1: Biot-DA 

The aim of this first test is to check the validity of the numerical method 
presented above using Biot-DA model. The domain is discretized with N x = 
700 which amounts to 32 points per slow wavelength and 104 points per fast 
wavelength, and iV = 6 diffusive variables are used. The source point emits 
symmetrically rightward and leftward moving fast and slow compressional 
waves, which are denoted Pf and P s , respectively, in figure EJ It can be seen 
from this figure that the numerical and analytical values of the pressure after 
200 time steps show excellent agreement. 
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-0.02 -0.01 0.01 0.02 0.03 2 3 4 5 6 7 

x(m) x(m) X 1Q- 3 



Figure 6: Test 1. Fast waves Pf and slow waves P s emitted by a source point at Xq = m. 
Comparison between numerical values (circle) and Biot-DA analytical values (solid line) 
of p at ti ~ 8.63 10 -6 s. Right: zoom on the slow wave. 



As mentioned above, the frequency correction introduced into the Biot- 
JKD and Biot-DA models has a negligible influence on the fast wave: con- 
sequently, we will therefore focus here on the slow wave. The error between 
the exact and numerical solutions will be measured in the L 2 norm in a do- 
main centered on the rightward-moving slow wave at time t\ ~ 8.63 10~ 6 s: 
[0, 0.007] m. Numerical values of the relative error and convergence order 
are summed up in table [3] at various values of N x and given in figure [3- (a). 
The convergence rate obtained by performing a linear regression is 2.00885, 
which is very similar to the theoretical second-order of the global algorithm. 

Figure 0-(b) shows the computational time in terms of the number of 
diffusive variables N, with N x = 700. The complexity of the scheme in term 
of diffusive variables is found to be in 0(N 2 ). 

With N = 6 diffusive variables, the linear optimization procedure de- 
scribed in section 13.31 yields: a\ = —588.77, a 2 = —365.69, a 3 = 369.78, 
04 = 1247.63, a 5 = —1956.56 and a 6 = 5725.45. Since some of the coef- 
ficients are negative, one cannot confirm that E is a decreasing energy in 
proposition |5J To examine this question numerically, the time evolution of 
E 3 in (13 3 p and —dE / dt in (I34p is shown in figure 8, where it can be seen 
that E 3 > 0, hence E > 0, and that dE / dt < 0. Despite the negativity of 
some ci£, figure 8 indicates that E is a decreasing energy and that Biot-DA 
is a well-posed problem. 
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N x 


Error L2 


Order 


1000 


2.23510" 


-1 




2000 


1.57110" 


-2 


3.853 


3000 


4.17610" 


-3 


3.353 


4000 


2.04710" 


-3 


2.558 


5000 


1.27310" 


-3 


2.166 


6000 


8.81710" 


-4 


2.032 


7000 


6.49410" 


-4 


1.989 


8000 


4.989 10" 


-4 


1.977 


9000 


3.95310" 


-4 


1.974 


10000 


3.21010- 


-4 


1.975 



Table 3: Test 1: error measurements and convergence orders. 




Figure 7: Test 1: relative error between exact and numerical solutions (left) in terms of 
the number of grid nodes N x . The dashed line is proportional to iV- 2 . CPU time (right) 
in terms of the number of diffusive variables N. 
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t(s) x1(f 5 t(s) x1(f 5 



Figure 8: Testl. Left: time-dependent evolution of E3 (|33|) : right: time-dependent evolu- 
tion of -dE I dt (J3U). 



5.3. Test 2: Biot-JKD 

The aim of the second test is to check the validity of the mathematical 
and numerical methods used to approximate the physical Biot-JKD model. 
Figure compares the numerical pressure obtained with the Biot-DA model 
with the analytical pressure obtained with the Biot-JKD model, at times t\ 
and t 2 > t\. The dispersion and attenuation of the slow wave can be clearly 
observed. Excellent agreement is found to exist between the two solutions. 

Two errors should be mentioned here: the modeling error e m , defined as 
the difference between the Biot-DA and Biot-JKD models; and the numerical 
error, e n , resulting from the numerical discretization of the Biot-DA model. 
The total error e t obviously satisfies: 

£t<£ m + £n- (49) 

Based on section [3731 taking N = 6 yields e m = 5.58 %. In test 1, e n ~ 1.70 
% was measured. At £1, the total error e t = 1.95 %, which means that the 
inequality f H9|) is satisfied but not optimally: the overall results are more 
accurate than those predicted on the basis of the bound f )49|) . The results of 
this test confirm that the method presented above efficiently approximates 
the transient waves modeled by the Biot-JKD model. 

5.4. Test 3: variable medium 

The aim of the third test is to establish whether the numerical methods 
presented in this paper can be used to handle more complex media. As an 
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x(m) 
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(d) 




0.007 0.008 0.009 0.01 
x (m) 



0.011 0.012 



Figure 9: Test 2. Fast waves Pf and slow waves P s emitted by a source point at t\ ~ 
8.63 10~ 6 s (a-b) and ti ~ 1.51 10~ 5 s (c-d). Comparison between the numerical Biot-DA 
pressure (circle) and the exact Biot-JKD pressure (solid line). Right row: zoom on the 
slow wave. 
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example, we took the porous medium with the parameters defined in table 
[21 except for the ratio tj/k, which varies linearly from 1.5 10 4 Pa.s.m -2 at 
x = —0.04 m to 5 10 9 Pa.s.m -2 at x — 0.04 m. These values are purely 
numerical and are not based on real data. Some changes had to be made to 
the method in comparison with that used in the homogeneous case: 

• at a given level of accuracy e m , the most-penalizing number of diffusive 
variables N has to be determined; 

• the coefficients ag have to be computed and stored at each grid point. 

In f )29p . the diffusive matrix S therefore differs between the grid points. In 
this example, the propagation matrix A remains unchanged since only the 
diffusive part is modified. When dealing with a real continuously variable 



medium, which occurs in the case of many applications [15[, the present 
ADER scheme would also have to be modified in order to handle the spatial 
changes in the matrix A accurately. 



250- P 




-0.02 -0.015 -0.01 -0.005 0.005 0.01 0.015 0.02 
x(m) 



Figure 10: Test 3. Pressure at ti ~ 8.63 10 6 s (zoom on the slow wave). 

Figure ffTUl) shows the pressure p at t\ ~ 8.63 10~ 6 s. As was to be 
expected, the rightward-moving slow wave is more strongly attenuated than 
the leftward-moving one, because the values of t]/k are higher in the right 
part of the domain. The present numerical tool therefore provides useful 
means for computing solutions of this kind, where no analytical expressions 
are available. 
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5.5. Test 4-' a 2-D example 

The one-dimensional method presented here can easily be extended to 
other dimensions. As a preliminary example, we take a two-dimensional 
medium with the parameters given in table [2j The number of physical un- 
knowns increase in this case from 4 to 8, and the equations of motion are 
also written in the form of a first-order hyperbolic linear system. The prop- 
agative part is solved with the ADER 4 numerical scheme. The diffusive 
part involves an order 1/2 fractional derivative for each component of the 
filtration velocity. The computational domain is set at [— 0.08, 0.08] 2 m. A 
Ricker source point, with a central frequency of 200 kHz and a time shift 
10~ 5 s, is localized at point (0, 0) and applied to the a xy component of the 
stress tensor. Applying our method with N = 6 diffusive variables to a grid 
of N x = N y = 1400 points gives the results presented in figure [TTJ Fast and 
slow compressional waves are observed as regards the pressure, while the 
additional shear wave is present in the a xx component of the stress tensor. 
It is proposed in future studies to address the analytical solution of the 2D 
Biot-JKD model and to perform an error analysis of the results obtained 
with the Biot-DA model. 





Figure 11: Test 4. Graph of the pressure (left) and the stress component a xx (right) 
emitted by a source point, at time t = 8.54 10~ 6 s. 



6. Conclusion 

A numerical method is presented here for simulating transient poroelastic 
waves in the high-frequency range. The Biot-JKD model, which involves or- 
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der 1/2 fractional derivatives, was replaced here by an approximate Biot-DA 
model, which is much more tractable numerically. Contrary to the approach 
used in [24), the Biot-DA coefficients are determined here using an opti- 
mization procedure, which depends on the frequency range of interest. The 
number of parameters and the accuracy of the model were quantified. The 
hyperbolic system of partial differential equations was discretized using effi- 
cient tools (Strang splitting and the fourth-order ADER scheme). The sta- 
bility condition of the numerical scheme is always independent of the param- 
eters involved in the approximate Biot-DA model. Numerical experiments 
performed in some academic cases (1-D homogeneous media) confirmed the 
reliability of this approach, and some preliminary simulations (with variable 
media, or in the 2-D context) show that the method is applicable to complex 
media. 

Some suggestions for future lines of research: 

• Heterogeneous porous media. Methods of modeling material interfaces 
in the context of Cartesian grids have been previously developed, based 



on an immersed interface method [22J. The possibility of applying this 



method to porous media in the low frequency range was studied in 
@, 0, 0, |27| ■ Work on means of extending this method to the Biot- JKD 
model is currently in progress. 

Thermic boundary-layer. In cases where the saturating fluid is a gas, 
thermo- mechanical effects have to be taken into account. Extended 



versions of the Biot- JKD have been developed [20j , involving additional 
order 1/2 fractional derivatives. The numerical method developed in 
this paper should lend itself well to working with this model. 
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Appendix A. Proof of proposition [JJ 

The equation (I9al) is multiplied by v s and integrated 
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The first term in (lA.ip is written 



dt dt 2 

Integrating by part and using (|9dl) . we obtain 

da f dv s 

v s — dx = / ——adx, 
Ox Jm ox 



dv s d 1 , 
pv s —-dx = — - \ pv s dx. (A. 2) 



/ ^7 (Ce - pp) dx, 
Jn dt 



de , f „ de 



C£ m dx -V p Tt dx ' (A3) 

jt\L Ce2ix -Sj p Tt ix - 



/3p^-dx. 
at 



The equation (|9b|) is multiplied by u> and integrated 



Pjw — + p w w 



^ + !L _L w (D + n)V 2 w + w^-) dx = 0. (A.4) 
dt k y/U dx) 



The second term in ( 1A.4I) can be written 



9* dt 2 

Integrating by part and using f )9d|) . we obtain 

<9p /" dw 

w — — dx = — / p — — dx, 
ox Jn ax 



p w w——dx = - — / p w w 2 dx. (A. 5) 



P Tt dx, 

p Tt& + fie ) dx > (A - 6) 

1 dp f de 

— p — dx+ / (3p — dx, 
m dt J R at 

— [ — f — jp 1 dx j ~\~ f j3 p — dx . 
dt \2 J R m J J K dt 
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After adding (lA.ip and the first term in (IA.4D . there remains 



dw 



dv. 



dx 



d 
dt 



Pf v s w dx. 



(A.7) 



Equations ( 1A.lj) -( tA4|) and the diffusive representation (!T5|) yield 



J t (E 1 + E 2 ) 




7] 1 



6m + K 7T y/Q e 



wip d9 dx. 



To calculate the right-hand side of (IA.8j) . equation ( 1T7|) is multiplied by w or 



dt 



dt 



(A.9a) 
(A.9b) 



After performing some algebraic operations on ( lA.9bl) . ( 1A.9ajl and ( 1A.8I) . one 
easily obtains the relation ( |T9|) . It remains to prove that E is a positive 
definite quadratic form. This is obviously so for E 2 and E 3 . Concerning E\, 
we write 



A 



1 2 1 2 
2P V s + ^PwW + PfV s W, 

1 rp 

-X T HX, 



(A.10) 



where 



x = 




H = 


f P 


Pf ) 




K w J 




V Pf 


Pw ) 



(A.ll) 



Taking S and V to denote the sum and the product of the eigenvalues of 
matrix H, we obtain 



V = det H = pp w - p) = X > 0, 
S = TrH = p + p w > 0. 



(A.12) 



The two eigenvalues of H are therefore positive, which proves that A is 
definite positive and completes the proof. 
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Appendix B. Proof of proposition [3] 

From ( 12tjj) . the system of diffusive evolution equations writes 



N 



dv s p f ^ 

-XT = — 7 y. at Wi 
at p f— ' 



9w 



iV 



-7 a £ ^ 



9a 

7*=°' 

^ = 
dt ' 



(B.la) 

(B.lb) 
(B.lc) 
(B.ld) 



dt 



N 



Qw -^^aeipe- (6j + Q)ifjj, j = l,...,N. (B.le) 



e=i 



Equation (IB. lbh is multiplied by w and (IB. lej) is multiplied by ipj 

N 

= — 'v in 

dt 

----- it ir r j — i' j 2_ 



dw \ ^ . 

w — — = —7 w y. a e Wt 



(B.2a) 



Summing (lB.2aj) and (lB.2b|) gives 



fl w i>j - 7 a/ Vfc - (fy + ^ J = 1, iV(B.2b) 



9w , dtpj 



N 



Slwif)j-(6j+ti)if)j-'y(w+if)j) y]cuil>i, j = l,...,N. 

1=1 

(B.3) 

The left-hand-side of (lB~3|) is equal to | ± (V + ^J). Then f lBTbj) is multi- 
plied by and (IB.leP is multiplied by w 



w 



N 



^3 = -T&j ^ae^Pe, j = l,...,N, (B.4a) 

1=1 

^ = Qw 2 - 7 w ^a<V/-(^ + fi)«'V'i. J = l,-,A(B.4b) 
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Summing (lB.4ap and (1B.4bl) gives 



w 



di/jj 



dt 



-A 



dw 
~dt 



N 



Qw 2 — (9j+Q) wipj — '-y (w + ipj) ai ipe 



1=1 



j = l,...,N. 
(B.5) 

The left-hand-side of (1B.5|) writes (w ipj). Elementary calculations on (1B.5j) 
and dEl yield (j = 1,...,N) 

d (\ 



dt \2 
Taking 



(w + if;, ) — w ip. 



(n w 2 - (6j + 2 n) w ipj + (6j + n) v>J) . 



E., 


1 , 
— - (w - 

2 v 

N 


E 


3=1 


X, 


j W 

= V 


Hi 


( n 



(w-^) 2 >0, 



20) 



9j + n, 



and summing the relations (IB. 61) for j 



N yields 



(B.6) 



(B.7) 



N 



3=1 



(B.8) 



Since the matrix H,, is triangular, its two eigenvalues are Q > and 6j + Q > 
0. The quadratic form XjHjXj is therefore definite and positive, which 
means that the left-hand-side of (1B.8j) is strictly negative. The energy E 



derived from system (IB. II) is therefore decreasing, and hence the system 
(IB. II) is well-posed. 
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